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The aim of this paper is to introduce a new Monte Carlo method 
based on importance sampling techniques for the simulation of stochas- 
tic differential equations. The main idea is to combine random walk 
on squares or rectangles methods with importance sampling tech- 
niques. 

The first interest of this approach is that the weights can be easily 
computed from the density of the one-dimensional Brownian motion. 
Compared to the Euler scheme this method allows one to obtain a 
more accurate approximation of diffusions when one has to consider 
complex boundary conditions. The method provides also an inter- 
esting alternative to performing variance reduction techniques and 
simulating rare events. 



1. Introduction. Monte Carlo methods are sometimes the unique alter- 
native used to solve numerically partially differential equations (PDE) in- 
volving an operator of the form 



The operator L is the infinitesimal generator associated with the solution of 
the stochastic differential equation (SDE) 



(1) X t = X + / a{X s )dB s + / b(X s )ds with era* = a. 
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It is well known that, for T > fixed, the solution on the cylinder [0,T] x D, 
of the parabolic PDE, 

( du(t,x) , 
-^-± + Lu{t,x) = Q, 

u(T, x) = g(x), for x € D, 

< u(t,x) = (f)(t,x), for (t, x) G [0, T]xdD, 

can be written as 

u(t, x) = E ttX [g(X T );T < r] + E^r, X T ); r < T], 

where r stands for the first exit time of X from the domain D. means 
that the process X is starting from x at time t. Thus an approximation of 
u(t,x) can be obtained by averaging g(XT)lT<r and 4>(t, X t )1 t< t over a 
large number of realizations of paths of X. Elliptic PDE may be considered 
as well. 

A large spectra of methods has been already proposed in order to simulate 
X (see, e.g., the books of Kloeden and Platen [22] and of Milstein and 
Tretyakov [29]). Most of these methods are extensions of the Euler scheme 
which provides a very efficient way to simulate (1) in the whole space. This 
method becomes harder to set up in a bounded domain, either with an 
absorbing or a reflecting boundary condition. Nevertheless some refinements 
have been proposed (see, e.g., [5, 15, 16, 19, 32, 34]). To improve the quality 
of the simulation or to speed it up, variance reduction techniques can be 
considered (see, e.g., [1-3, 17, 20, 21, 31, 38]). This list is not intended to 
be exhaustive. 

In the simplest situation, for a = Id and b = 0, the underlying diffusion 
process is the Brownian motion. Muller proposed in 1956 a very simple 
scheme to solve a Dirichlet boundary value problem. This method is called 
the random walk on spheres method [30]. The idea is to simulate successively, 
for the Brownian motion, the first exit position from the largest sphere 
included in the domain and centered in the starting point. This exit position 
becomes the new starting point, and the procedure is iterated until the exit 
point is close enough to the boundary. Nevertheless, simulating the exit time 
from a sphere is numerically costly. In [27], Milstein and Rybkina proposed to 
use this scheme for solving (1) by freezing locally the value of the coefficients. 
In a first approach, spheres (that become ellipsoids) were used. Later on 
[26] (see also [29]), Milstein and Tretyakov used time-space parallelepipeds 
with a cubic space basis. For this last approach, it is easier to keep track 
of the time but the involved random variables are costly to simulate. In 
order to overcome these difficulties, one may think to use tabulated values. 
This is memory consuming as the random variables to simulate depend on 
one or two parameters. The method of random walk on squares was also 
independently developed in the Ph.D. thesis of Faure [11]. For the Brownian 
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motion, this method is still a good alternative to the random walk on spheres 
(see [7] for an application in geophysics). 

In [8], we have proposed a scheme for simulating the exact exit time and 
position from a rectangle for the Brownian motion starting from any point 
inside this rectangle. Compared to the random walk on spheres method, this 
method has the following advantages: 

• It can be used whatever the dimension and, as for the random walk on 
squares, a constant drift term may be added. 

• The rectangles can be chosen prior to any simulation, and not dynamically. 
There is no need to consider smaller and smaller spheres or squares when 
the particle is near the boundary. 

• The method can be also adapted and used for the simulation of diffusion 
processes killed on some part of the boundary. 

The method we propose here is based on the idea to simulate the first 
exit time and position from a parallelepiped by using an importance sam- 
pling technique (see, e.g., [12, 14]). The exit time and position from a par- 
allelepiped for a Brownian motion with locally frozen coefficients is chosen 
arbitrarily, and a weight is computed at each simulation. By repeating this 
procedure, we get the density on the boundary or at a given time of the 
particles, by weighting the simulated paths. As we will see, the weights are 
rather easily deduced from the density of the one-dimensional Brownian mo- 
tion killed when it exits from [—1,1]. All involved expressions are numerically 
easy to implement. 

This new algorithm is slower than the Euler scheme for smooth coeffi- 
cients, but it is faster than the random walk on squares [7, 29] and the 
random walk on rectangles [8]. It can be used to simulate the Brownian 
motion as well as solutions of stochastic differential equations for specific 
complex situations as: (a) complex geometries (the boundary conditions are 
correctly taken into account); (b) fast estimation of the exit time of a do- 
main for the Brownian motion (only few rectangles are needed); (c) variance 
reduction; (d) simulation of rare events. 

This algorithm could be relevant for many domains: finance, physics, bi- 
ology, geophysics, etc. It may also be used locally (e.g., it can be mixed with 
the Euler scheme and used when the particle is close to the boundary) or 
combined with other algorithms, such as population Monte Carlo methods 
(see Section 4.5). 

We conclude this article with numerical simulations illustrating various 
examples. It has to be noted that choosing "good" distributions for the exit 
time and position from a rectangle is not an easy task in order to reduce the 
variance. We then plan to study in the future how to construct algorithms 
that minimize the variance, as in [1, 3]. We have to consider for this a high- 
dimension optimization problem. 
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Outline. In Section 2, we present the importance sampling technique 
applied to the exit time and position for a (drifted) Brownian motion from 
a rectangle. In Section 3, we recall briefly some results about the density of 
the one-dimensional Brownian motion with different boundary conditions. 
The explicit expressions are given in the Appendix. In Section 4, we present 
our algorithm and compute its weak error. Four test cases are presented in 
Section 5. We compare also our algorithm with other methods in this last 
section. 

2. Algorithm for the exit time and position from a right time space par- 
allelepiped by using an importance sampling method. The aim of this part 
is to give a clear presentation of our method. In order to avoid ambiguous 
notation we consider in this section the situation of a two-dimensional space 
domain. The results can be easily generalized to higher space dimension. 

We are looking for an accurate approximation of the exit time and position 
from a right time-space parallelepiped which is a geometric figure in the 
three-dimensional space. 

For L\,L 2 > given let R be the rectangle [—Li,Li] x \—L2,L 2 \. The 
rectangle R is the space basis of the right time-space parallelepiped Rt = 
[0, T] x R for a fixed T > 0. We can also consider R^ = M + x R, and set in 
this case T = +oo. 

For T < +oo, the right time-space parallelepiped Rt has six sides which 
are denoted by 



In other words, each side of Rt is labeled by a couple (i,rj) G {0,1,2} x 
{—1,1}. For i G {1,2} the side Si >r) is perpendicular to the unit vector in 
the ith direction. For i = 0, the side So,-l corresponds to the rectangular 
initial basis while the side So,i corresponds to the top of the time-space 
parallelepiped Rt for T < +oo (see Figure 1). 

From now on, we shall identify each side with the corresponding (i,r])- 
indices. 

We consider a time-homogeneous diffusion process (X t )t>o living in R. 
On each side of R, the process X may be reflected or absorbed. Moreover, if 
T < +oo, the process is stopped at time T. We can thus identify the sides of 
R with the sides Si jT) of Rt for i G {1, 2} and rj G {—1, 1}. We denote by SH 
the subset of {1,2} x {—1,1} that contains the indices of the sides on which 
a Neumann boundary condition holds (possibly, 9t = 0). On this set the 




m x r, 

{0} x R, 

[0,T] x [-Li,Li] x{ v L 2 } 
[0,T] x{7?Li}x [-L 2 ,L 2 ] 



for j) G { — 1, 1}, 
for 7] G {-1, 1}. 
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diffusion is reflected. Let us denote by D the subset of {1, 2} x {—1, 1} that 
contains the indices of the sides on which a Dirichlet boundary condition 
holds. On this set the diffusion is killed. Finally let us set 21 = 3D if T = +oo 
and 21 = {(0, 1)} U D if T < +oo. With this notation the time-space process 
1 1 — y (t,Xt) is absorbed when hitting one of the sides Si ;V with (i,rj) G 21. 

Let B = (B 1 , .B 2 ) be a two-dimensional Brownian motion and fi = (/ji, 112) 
a vector of M 2 . For % G {1, 2}, we set 

J 1, if (i,rj) G9t (reflection), 
^ ~ \ 0, if (t, 77) € 01 (absorption). 

We consider the two-dimensional diffusion process (X,¥ x ) xe ji whose co- 
ordinates are, for x = (xi,Xz) € R, 

(2) Xt = x i + Bl + i H t + li ^{X i )- li ^ 1 iJ Li {X i ), P x -a.s., 

where lf Ll (X l ) stands for the symmetric local time of X 1 at ±Lj, respec- 
tively. 

We define r = T,n = inf{< > 0| > for t G {1,2} and 

r = min Tj. 

i€{0,l,2} 

In addition, we set J = argminj g | 0i i,2} T %- With this notation, unless J 7^ 0, 
the Jth component of X is the first to exit from the domain. For J G {1,2}, 
let us define e = X^/Lj G { — 1, 1}. For J = we set e = 1. In this case X 
has not reached the sides of D before time T. 

The couple (J,s) labels the side in 21 of the parallelepiped Rt = [0, T] x R 
that the diffusion X hits first. Note that with our convention, the sides on 




Side So,i Side S 1 ^! Side Si t ^i 




Side So,-i Side S2. 1 SideS^-i 



Fig. 1. Convention for the sides of Rt = [0,T] x R. 
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which the process is reflected cannot be reached so that Tj = +00 if X 1 is 
reflected both at —Li and Lj. 

We are interested in computing E x [/(r, X T )\ by a Monte Carlo method 
for a bounded, measurable function / where r is defined as above. 

Instead of simulating (t,X t ), we will simulate some random variables 
according to the following procedure. The aim is to simulate (J,e,t,X t ) by 
using an importance sampling technique. In order to do this we choose a 
probability F x which is absolutely continuous with respect to ¥ x , and we 
draw a realization of (J,e,t,X t ). Let us set 

a itV =F x [(J,e) = (i,r))\ 

for (i, rj) G 21. For (i, rj) G 21 let ki^ denote the density under F x of (t,X t ) 
given {(t,X t ) G S ii?? }. 

In order to simplify notation let us consider an underlying probability 
space (Q,J-,F X ) rich enough. Let Z be a random variable on this space, 
with distribution F x . Let A be a measurable event on this space. We suppose 
that, conditionally on A, Z has a density with respect to the Lebesgue 

measure. Let us introduce the following convention: 

F x [Z = z;A]=p(z\A)F x [A]. 

That is, for B a measurable event of (ft, J-, ¥ x ), 

F x [{Z€B}nA}= ( p{z\A)¥ x [A]dz= ( F x [Z = z;A]dz. 

J B JB 

Consider now the following notation: let (z, 77) G 21. For i G {1, 2} set j = 3 — i. 
Then for any 6 > and z G Si^, we define 



(3) M itV (e,z) 



x [n = 9; X* T . = ijL^ [X 3 e = Zj ; Tj > 
ati >v ki,ri(Q,z) 



where kij, is the {X T G S'j^j-conditional density under F x of (t,X t ). 
If T < +00, we define 

(4) M 0)1 (T,z) = 1 J] Vx[xi. = z j ; Tj >T\, 

where is the {X T G Sj^j-conditional density under F x of (r, X T ). 
We call Mi jT} weight. 

Proposition 1. The weights M i)V defined in (3) and (4) satisfy 
E x [f(r, X T )} = E x [Mj,(t, X T )f(r, X T )] 
for any measurable and bounded function f on 8Rt ■ 
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Before proving this proposition let us introduce the algorithm. 
The algorithm is described as follows: 

Algorithm 1. Let x be fixed in R. 

(1) Draw a realization (J,e) of (J, e) E 21 under F x . 

(2) Draw a realization of the exit time and exit position (r, X-) according 
to the density k-j- on Sj-. 

(3) Compute the value of Mj-(r, X T ) by 

E x [M J>E (r, X T )f( T , X T )\ = E x [f(r, X T )\ . 
We call Mj-(t,X t ), weight. 

If {(J ,eW,rW,X~ j^W)}^!,...^ are iV independent realizations of the 
random variables ( J, e, r, X T , Mj t£ (r, X T )) constructed as above, by the law 
of large numbers we have 

1 N — 

E x [f(r, X T )] = - ^ «J« /(r« , X% ) . 

°° i=i 

The main feature of our approach is that the weights Mj j£ (r, X T ) can be 
easily evaluated. 

Remark 1. In order to evaluate with (3) and (4), there is no need 
to know ¥ x [(J,e) = (i,rj)]. It is important to notice that Mi v depends only 
on the one-dimensional distributions of the drifted Brownian motion. 

Proof of the Proposition 1. We want to prove that 

E x [f(r, X T )\ = E x [Mj, £ (t, X T )f(r, X T )} 

for any measurable and bounded function / on ORt- 

We remark first that if pi tV = ¥ X [(J, e) = (i,r/)] for (i,rj) in 21, then 

(5) E x [f(r,X T )}= ^ L %[M i , v (r,X T )f(T,X T )\(J,e) = (i,r,)]. 

Furthermore, for (i, rj) £ 53, if i = 2 set j = 1 and z = (zi,r)L2) else, if i = 1 
set j = 2 and z = (r/Li, 22)- 

E a: [/(r,X T )|(J,e) = (i,r ? )] 

= / /(0 )Z )P a: [(Ti,X4) = (0,^)l(^e) = ^'7)]^^ ) 

J[0,r|x[-£ J -,Li] 
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where P x [(r,, X 3 Ti ) = (9,Zj)\(J,e) = (i,r])] is the {(J, e) = (i, ^-conditional 
density of (ri,Xr t ) with respect to dtdzj. Hence 



E x [f(r,X T )\(J,e) = (i, V )]=E x [f(r,X T )Ml v (T,X T )\(J,s) = (i, V )], 



where 



MUe,z) 



Let us note that Mi iV (9,z) 
that 



[(r i ,X^) = (e,z J )\(J,e) = (i, v )] ^ 

^i,r)( T 5 X T ) 

- M[ (9,z)pi >r) /ai ir) . With (5), we can deduce 



E x [f(r,X T )]=E x [f(r,X T )M J>e (T,X T )]. 
Indeed, it suffices to remark that for (i,rf) G J), 



M^(6,z) 



ai jV k ijri (9,z) 



1 



: [(r i ,X^ = (e,z j );(J,s) = (i,r,)] 
,[{Ti,X^) = (0,z j );Xi i =r,L i ,r^> 



a iin k itV (6,z) 

The independence of the coordinates of X leads to the desired equality. If 
T < +00, similar computations imply that for z G R, 

1 



M 0A (T,z) 



ao,ik ,i{T,z) 
and the conclusion also holds. □ 



Xt = z; min n > T 

iG{l,2} 



(6) 



Let us evaluate these probabilities. 
For i G {1,2}, let p l (t,x\,x 2 ) be the solution of 

f dp l (t,xi,x 2 ) 1 d 2 p l . . dp 1 

it Ti To I -A- II; 

dx 2 

for (t,xi,x 2 ) el+x (-Li, Li) 2 , 



Of 



2 dx 2(^ x ^ x 2) + l^i— (t,xi,x 2 ) 



p t {t,x 1 ,x 2 )-^5 Xl (x 2 ), for (xi,x 2 ) G (-Li, Li) 2 , 



with the following boundary conditions (b.c): 

p i (t,x 1 ,-L i ) = (Dirichlet b.c.) if (*,— 1) G 2t, 
dp 1 

— — (t, x\, —Li) = (Neumann b.c.) if (i, —1) G D\, 

ox 2 

p i (t,xi,L i ) = (Dirichlet b.c.) if (i, 1) G 21, 

dp 1 

— — (t, x\,Li) = (Neumann b.c.) if (i, 1) G 1H. 

ox 2 
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Thus, p % denotes the density of the drifted Brownian motion X 1 with possibly 
some reflection at the endpoints of {—Li, Li), and killed when it exits from 
this interval by an endpoint where no reflection holds. For / a bounded 
measurable function from [—Li, Li] to R, we have 

/Li 
p\t,x l ,x 2 )f(x 2 )dx 2 
-Li 

for x\ £ [—Li, Li] where P Xl is the distribution of X % with Xq = x\ € [—Li, Lj\. 
Let us note that the distribution of the marginal X 1 of X under P( xliX2 ) de- 
pends only on x%. 

We introduce the scale function <3> J ' + of X 1 defined by 



^/2[LiLi — g l^liLi ' ' 
X2 + Li 

if m = 0. 



for x 2 € [-Li, Li], $ i '+( X2 ) = < 

~2Li ' 

The function & l ' + (x 2 ) has been normalized such that $ l,+ (Lj) = 1. Let us 
note that & l,+ (xi) = P^fX* = Li] if Dirichlet boundary conditions hold at 
both endpoints — L, and Lj. We also set & % ~{x 2 ) = 1 — <& l ' + (x 2 ). 

If Dirichlet boundary conditions hold both at —Li and Li, then we set for 
t > and (x\,x 2 ) £ [-Li, Li] 2 , 

$*' ± (x 2 ) 

p 1 ' (t,X 1 ,X 2 ) = p\t,X 1 ,X 2 ) ± . 

Via a Doob transform, for a bounded and measurable function /, 

E Xl [f(Xi);t < r,|X; = ±Li] = p i > ± (t,x 1 ,x 2 )f(x 2 )dx 2 . 

J -Li 

Let us set for x\ G (—Li, Li), 

f Ll dp i 

(7) q\t,xi) = - —(t,xi,x 2 )f{x 2 )dx 2 

J —Li 



q i > ± (t,x 1 )=- I ^(t, Xl ,x 2 )f(x 2 )dx 2 . 



and 

We can easily deduce that 

Pa* fa < i] = / g*(a, xx) ds and P^ [t< < t\X* n = ±L t ] = [ <f '±1 

JO ■/ 



s, x±) ds. 



In other words, q l (t,x\) [respectively, q l,± (t, x\)] is the density of the first 
exit time from [—Li, Li] for X 1 (respectively, the first exit time from [—Li, Li] 
for X i given {X\. =±Li}). 
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Thanks to these expressions, Mq^(T, z) and Mi tTI (9, z) are easily computed 
since 

p\e,xi,zi), 

q i ' ± {9, x l )¥' ± {x i ) if (i, -1) G 21 and (t, 1) G 21, 
g^Mi) if (i, —1) G £H and (i, 1) G 21, 
g^Mi) if (i, — 1) G 21 and (i, 1) G £H. 

3. Analytical expressions for the densities. In order to compute p % (t, x\, X2), 
together with q l (t,x\) and q 1 ' (t,x±) by (7) and (8), one has to solve equa- 
tion (6). By using a scaling principle, we may assume that Lj = 1, as 

i ( , , I f t xi x 2 \ 

p (t,xi,x 2 ) = —pijz, —, —; Li/j, j, 

where p(t,x\,X2',5) is solution to (6) with Lj = 1 and a convective term m 
equal to 5. 

There are basically two ways to obtain p(t, x±,X2',5). The first one is based 
on the spectral expansion of iA + <5V since this operator may be reduced to 
a self-adjoint one with respect to the scalar product induced by the measure 
exp(— 25xi). The second one is the method of images when 5 = 0. 

If 5 7^ 0, the case of a Dirichlet boundary condition at both endpoints may 
be treated by using a simple transform that reduces the problem to 5 = 0. 

For the case of Neumann boundary condition at both endpoints, one can 
invert term by term the Laplace transform of a series for the Green function. 

In the case of a mixed boundary condition, the previous method gives 
rise to a series that cannot be used in practice, so only the spectral expan- 
sion should be used. In addition, the first eigenvalues have to be computed 
numerically. 

As the formula are standard in most of the cases, we give the relevant 
expressions in the Appendix. 

4. General domain. As stated before, we aim to solve by a Monte Carlo 
method a parabolic or an elliptic PDE. The idea is to represent the domain 
as the union of time-space parallelepipeds and to simulate the successive 
exit times and positions from these parallelepipeds. Attention has to be 
paid while doing this decomposition in order to control the error at each 
simulation step. 

4.1. From parallelepipeds to right parallelepipeds. Consider herein the 
notation of Section 2. Let us study first the parabolic PDE with constant 



P Xi [X % e = z i ]T i >T\ = 
f> Xi [r i = e-Xl = ±L i ] = 

F Xi [T l = 6;X l e = L i ] = 
F x .[n = 9;Xi = -L i ] = 
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coefficients A, c and /i = (/ii)i=i,...,d on the rectangle Rt, 

dv(t,x) l-s^d 2 v(t,x) 
+ 9 Z-^i ' 



II 



(9) 



dt ' 2 

d 

+Z> 



x) 



i=l 

dv(t, x 
dxi 



dxf 



+ cv(t,x) = \, on Rt, 



0. 



for x E Si jJ? if (i, r/) E 



for x E Si jT) if (i, r/) E 21, 
if T < +oo. 



<9xj 

u(i,x) = cf>(t,x), 

We assume that a classical solution to this problem exists, which is, for 
example, the case if <f> and g are continuous and bounded. Let X be the 
diffusion process whose components are given by (2). Then it follows from 
the ltd formula applied to X that, for t E [0,T], 

v(t, x) = E x [e c ( r -'V(r - t, X T _ t ); r < T - t] 

rr-t 



+ E x [e c ( T - t ) 5 (X T _ i ); r = T-t]+E a 



X 



e c(r-t-s) ds 



where r is as above the first exit time from Rt- 

Let us remark that if a is an invertible d x ti-matrix, then the function 
u(t,x) = v(t,a~ 1 x) is solution to 

du(t,x) i 1 ^ r __*i d 2 u(t,x) 
2 



(10) 



(7 <J 



a 



i=l 

chx(i, x) 



<9xj 



0. 



J ' 1 5xj 
■u(i, x) = cr x x), 
_ «(T,x) =g(cr- 1 x), 



on [0, T] x 

for x E if (i, r/) E 

for x E crSi :rj if (z, 77) E 21, 
if T < +oo. 



If n.j is the unit vector orthogonal to the side crSi^, then rij = (o*)~ l e; L , 
where is the unit vector in the ith direction. It follows that <7<7*rij = ere,/ 
and thus 

r i- q \ *\ xj u \ du(t,x) 
for x E abi ±\ [era rij • Vu{t, x) = cr ? - j — , 

OXj 

which means that a Neumann boundary condition in the co-normal direction 
holds in (10) on aSi jT] if (i,rj) E 91. 

We can thus solve (10) by reducing the problem to (9) and use a Monte 
Carlo method in order to compute the values of u(t,x). 
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4.2. The hypotheses. Let us consider a domain Q in 



. For the 



sake of simplicity, we assume that Q is the cylinder [0, T]x D (with possibly 
T = +oo), where D is an open, bounded domain of R d with piecewise smooth 
boundary. Let us consider a function a with values in the space of d x 
(i-symmetric matrices which is continuous on D and everywhere positive 
definite, together with some functions b:Q^R d , c:Q^R and f:Q^R. 
For all (t,x) € Q, we denote by o~(t,x) a d x (i-symmetric matrix such that 
a(t,x)a*(t,x) = a(t,x). 
We set 



L 



1 d 

2 Yl a M 



(t,x) 



3 2 3 
+ ^bi{t,x) 

3 i=l 



dx-i dx 



dxi 



Let us introduce the hypotheses needed to ensure the convergence of our 
algorithm. To set up a Monte Carlo numerical scheme, one needs three inter- 
connected ingredients: 

(i) The existence and the uniqueness of a solution u to the following 
PDE 



' du(t,x) 
dt 



(11) 



+ Lu(t, x) 



+ c(t,x)u(t,x) + f(t,x) 
u(T,x)=g(x), 
u(t, x) = 4>(t, x), 
k d n u(t,x) = 0, 



on [0,T] x D, 
x G D, 

on T d C [0,T) x 3D, 
on T n C [0,T) x 3D, 



where d n denotes the co-normal derivative along the lateral surface. Td (re- 
spectively, T n ) are subsets of [0,T) x 3D on which a Dirichlet (respectively, 
Neumann) boundary condition holds. 

(ii) The existence of a solution to the diffusion process associated with 
L. Note that since the simulation involves distributions and not stochastic 
integrals, we do not need strong existence for the associated SDE. 

(hi) The solution u can be expressed in terms of the probabilistic repre- 
sentation 



u(t,x) = E 



t,x 



(12) 



+ E 



+ E 



t.x 



c(s, X s )dsj (f)(r, X T )t T<T 
exp^ c(s,X s )ds^jg(X T )t T>T 
expfy c(r,X r )dr)f(s,X, 



ds 



where r is the first exit time from [0, +oo) x D by a point of T^. 
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Notation 1. We denote by V the set of time-space parallelepipeds P 
such that there exist < s < t < T, L±, . . . , L d and x G M. d such that 

P = [s, t]x(x + a([-L u LJ x • • • x [-L d , L d ])), 

where a is & d x <i-matrix. Possibly t = +oo (if T = +oo). 

The assumptions that have to be done are the following: 

(HI) There exists a subset Vd of V such that Q = Upe-p D P- Besides, if P = 
[s, t] x U G V for a parallelepiped U, then for all r G [s, t), [r, i] x [/ £ P. 
In other words, one can truncate the parallelepipeds in time. 

(H2) There exist T n , contained in dQ = [0, T] x dD and some subsets 

Pa-, Pd of / such that T n C Upe^ ^^°) C UpePd ^ ne c l° sure of 
T n U is equal to [0,T] x 5L> and r n n = 0. This means that the 
boundary of [0, T] x dD is split in two distinct parts, where either the 
Dirichlet or the Neumann boundary conditions hold. More precisely 
a side of a parallelepiped in Vd contained in dQ is either from T n or 
from r,j. 

(H3) The differential operator L is the generator of a continuous diffusion 

process X that is reflected at T n and killed when hitting TdU {T} x D. 

The probabilistic representation of the solution given by (12) holds 

(see, e.g., [25] for existence results of such reflected process and [36] if 

there are no reflections). 
(H4) There exists an unique solution u of class C 1 ' 2 on [0,T) x D to (11) 

which is continuous on [0, T] x D. 
(H5) For a right parallelepiped R and a matrix a let P = [s, t] x {x + aR) G 

Vd- We associate with P a vector b G M. d , two constants c, / and we 

construct the differential operator 



-_i y ~ ki _?_ + Yb k — 

2 4^ ' dxudxf cL. 
fcj=i K * fe=i Xfc 



with a = era* . 

Fix 5 > 0. We assume that the solution u to (11) satisfies, for any 
y in the interior of x + aR, 



where X is the diffusion process generated by L, and f is its first exit 
time from P. 

Remark 2. If T = +oo and the coefficients are time-homogeneous and 
Td = [0, oo) x 7^, r n = [0, oo) x 7 n , then v(x) = u(0,x) is solution to the 
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elliptic PDE 

(Lv{x) + c{x)v{x) = f{x), on D, 

v(x) = 4>(x), on 7 d C 9D, 

d n v(x) = 0, on 7 n C 3D. 

Thus, by solving the parabolic PDE (11), we may also solve the elliptic PDE 
(13). We will thus focus only on (11). 

Remark 3. The result of the existence of a stochastic process reflected 
on some part of the boundary of [0, T) x D is deduced from the existence 
of a stochastic process reflected on the lateral boundary [0, T) x D which is 
killed when it hits T n . 

4.3. The algorithm and its weak error. In order to simplify the notation, 
if T < +oo, we denote the final condition g of (11) by <j){T,x). 

Given (t, x) € Q, the solution u{t, x) of (11) is computed by the Feynman- 
Kac formula. For this, we have to simulate the diffusion process X up to its 
first exit time r from Q. We suppose here that the particle cannot exit by 
a part of boundary where a Neumann boundary condition holds. Let u be 
the solution of (11). Let us introduce the following notation: 



for s > t 



Y s = 1 + jf c(r, X r )Y r dr = exp Qf c(r, X r ) dr^j , 
Z s = J" f(r,X r )Y r dr. 



Then u(t,x) is given by 

(14) u(t,x)=R tiX [<f>(T,X T )Y T + Z T ]. 

We construct now the algorithm that approximates (14) by a Monte Carlo 
method. 



Algorithm 2. Assume that we start initially at the point (t,x) G Q 
and fix a number N of particles. 

(1) For i = 1,...,N do 

(A) Set (0„, Ho, Y Q ,Z ,W ) = (t,x, 1,0,1) and k = 0. 

(B) Repeat: 

(a) Choose an element G Vd of the form P^> = [6k, s] x U, 
U C M d such that {6k, Hfc) belongs to the basis of P (s is pos- 
sibly infinite if, for example, T = +oo and the coefficients are 
time-inhomogeneous). On P^ k \ consider the differential opera- 
tor as well and which approximate L, c and / as 
in (H5). 
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(b) Draw a realization of a random variable (0k+i, ^k+i) with val- 
ues in ({s} x U) U ((9k, s) x dU) and compute its associated 
weight w k as shown in Sections 2 and 4.1 by considering the 
exit time and position from the parallelepiped P^ . 

(c) Compute W k = W k -\w k and 

y fe+1 = y fc exp( c ( fc )(# fc+1 -# fe )), 

Z k+ i = Z k + /(*> f k+1 exp(c( fc ) S ) ds. 

(d) If E k+ \ G Td or = T, then exit from the loop. 

(e) Increase k. 

(C) Set (9^,E^,Y^,Z^,W^) = (9 k+1 ,E k+1 ,Y k+1 ,Z k+1 ,W k ). 
(2) Return 



1 N 

(15) u(t,x) = — ^2(W {i U(0 {i \E^)Y^ + W®Z®). 

i=l 

We denote from now on by F x the distribution of the Markov chain A k = 
(9 k ,E k ),k > 0. Note that (Y k , Z k , w k ) k > is obtained from (A k ) k > . 

Proposition 2. For any (t,x) <E [0, T) x D, 

(16) \u(t,x) - E x [u(t,x)]\ < 6E x [W 1/ uexp(M9 u % 

where 5 is defined in (H5), v is the number of steps that (A k ) k >o takes to 
reach the boundary n {T} x D and 

M= sup c(s,y). 

(s,y)e[t,T)xD 

Remark 4. Note that the weak error in (16) does not depend on the 
choice of the importance sampling technique while the Monte Carlo error 
depends on this choice. If the coefficients a, b, f and c are constant on the 
domain, one can choose 5 = and the simulation becomes exact. 

Proof. To the Markov chain (A k ) k >o is associated a random sequence 
of parallelepipeds (P^) k =o,...,v Let us denote by the successive times 
the diffusion process X reaches the boundary of the P^ 's. 

Since Zq = 0, Yq = 1 and u = <j) on the boundary of Q, we get 

%[u(t,x)]=%\W v Y v <l>(e v ,3 v ) + W v Z v ] 



(17) =u(t,x) + E 2 



u-l 



W v y~]( z k+i ~ z k + Y k+1 u(9 k+ x,E k+ x) 



fc=0 



16 



M. DEACONU AND A. LEJAY 



Y k u(6 k ,E k )) 



Let (Qk)k>o be the filtration generated by the Markov chain (A k )k>o- 
We remark that Y k and Z k are measurable with respect to g k while w k is 
measurable with respect to Gk+i (since it is obtained from 6 k , E k , k+ \ and 
Sfc+i). 

By using the Markov property, after setting W k+ i )V = E x [w k+ i • • • w v \Gk+i], 
we get 

E x [W u (Z k+l -Z k )} 

= E x [W k+1)U E x [w k (Z k+1 - Z k )\g k ]W k -!], 
E x [W u {Y k+ iu(6 k+1 ,Z k+1 ) - Y k u(6 k ,E k ))) 

= E x [^ fc+li A[^(n+i«(^+i,s fc+1 )-y fcU (^,H fe ))|g fc ]w fc _ 1 ]. 

Let us denote by (X( fc ),pfj) the process generated by the operator 

with constant coefficients a" and on p( k \ Define recursively (t^°\x^) = 

(t,x) and (t (fc+1 ),x (fc+1) ) = (^,1"^],) where r (fc ) is, as above, the first exit 

time from PW for the diffusion X^ k \ Let also and be the values 
that approach / and c on P^ k \ and define also recursively = 1 and 

y {k) = y(k-i) exp( c (*)(t(* +1 ) - t( fc ))). 

By using the properties of ¥ x and the ltd formula we obtain 
E x {w k (Y k+1 u(6 k+1 ,E k+1 ) -Y k u(6 k ,E k ))\g k ] 



y 



(*) E (fc) 



04 



Also, 



E x [w k (Z k+1 -Z k )\g k ] = y^E { f k) 



f (k) 



t(*+i) 



„(*), 



Under the hypothesis on the coefficients and the parallelepiped P( fc ) we have 
|Ea;[^fc(^fc+i«(0fc+i,S fc+ i) - Y k u(O k ,E k ) + Z k+1 - z k )\g k ]\ 
r „t( fc+1 ) 



=c (*)( s -t(fc)) 



t(fc) 



at 



(k) 



u(s,lf')|/W j 



<y(%<E x [<^ fe y fe |g fe ], 
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since the Y^'s (and so the y^'s) are positive. Hence, from (17) and the 
Jensen inequality applied to | • | , we obtain 



\E x [Y u <p(6 u ,E u ) + Z v ] -E x [u(t,x)]\ < SE X ^ 



u-l 



As < Yfc < e AWk for k = 0, . . . , v, we deduce (16). □ 

4.4. The Monte Carlo error. In order to compute the solution u(t,x) of 
(11), we have constructed the estimator u(t,x) given by (15) whose variance 
is 

Var^ u(t, x) = jj VargjW^^, E V )Y U + W V Z V ). 

The Monte Carlo error depends on this variance s 2 = Var^ u(t,x), since 

asymptotically for iV — > oo the true mean E 2 ,.[ : u(t, x)] lies in the interval 
[u(t,x) — 2s,u(t,x) + 2s] with a confidence of 95.4%. 

We denote by P n the distribution of (Afc)^>o with respect to the real 
distribution of the exit time and position of the rectangles. In this case the 
weights are equal to 1. Any event $ measurable with respect to (Ak)k>o 
satisfies P n [$] = P X [W$]. 

We get thus 

V&i fx (W u <p(6 u ,Z u )Y u + W U Z U ) = * + Vax fn (u(t,x)) 

with 

* = E n [(W„ - l)((/>(e u ,E v )Y v + Z v )\ 

This shows that a good choice for the density of the exit time and position 
from the parallelepipeds is such that < is as small as possible. By the 
way, reducing the variance is a difficult task and requires some automatic 
selection/optimization techniques, as explained in the Introduction. 

In addition, the numerical experiments we performed up to now highlight 
another difficulty. W u may take large values, and this implies meaningless 
values for u(t,x). That is why we suggest to keep track also of the empirical 
distribution, or at least of the variance of W v . 

In order to illustrate this, let us assume that the diffusion process X has 
no drift and that for the simulation, the right parallelepipeds we use are 
squares centered on the particle, and consider the same density for the exit 
time and position. By a scaling argument, the distribution of the weight Wk 
at the feth step does not depend on the size of the squares, so that the w^s 
are independent and identically distributed under ¥ x . 

Let us fix an integer n such that v>n a.s. (for example, the minimal num- 
ber of steps needed to reach an absorbing boundary). We set £ l = log(u>j), 
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so that W n = exp(^™ =1 £*). As the £ l are independent and identically dis- 
tributed, let us note S n = Yli=i£ l > then S n /^/n converges to some normal 
random variable x with mean m and variance s 2 . For n large enough, the 
distribution of W n is close to the distribution of exp(- v /ra\ / )- We obtain, 
with the expression of the Laplace transform for the normal distribution, 
for je {1,2}, 

( fs 2 

^x[(W n y] KiE x [exp(js/nx)] = explmjy/n + n—^- 
This leads us to the following approximation: 

Varjj (W n ) w eyLp{2m\/n + 2ns 2 ) — exp (^rriy/n + — s 2 ^j 

m i\( ( m -A f m 3n 2 \\ 

« exp(2n, ) ^exp (_ + 1 j - exp - -s j j 



2\ 



~ exp(l + 2ns 

Tl— >00 

So, for large re, the variance of W n explodes, while E x [W n ] = 1 for any re > 1. 

In [13] (see also [14]), Glynn and Iglehart exhibit another argument that 
shows that the simulation performs badly if too many steps are used. 

4.5. Population Monte Carlo. In order to overcome the explosion of the 
variance due to the weights one can use a population Monte Carlo method. 
This kind of method, also known as quantum Monte Carlo, sequential Monte 
Carlo, Green Monte Carlo, . . . has been used for a long time in physical 
simulations (see, e.g., [18] for a brief survey) but also in signal theory, 
statistics, .... A probabilistic point of view is developed in the book [9] of 
Del Moral. 

In our case, instead of simulating the particles one after another, the 
idea is to keep track of the whole population of N particles (y^')ie{i,...,N} 
with time and space coordinates (t^\x^) and a weight u;W according to the 
algorithm given below. Each particle has two possible states: still running or 
stopped. A particle is stopped either at the first time it reaches an absorbing 
boundary, or if its time is equal to the finite final time T. Otherwise, the 
particle is still running. 

Algorithm 3. This algorithm computes an approximation of the quan- 
tity E x [f(T A r, Xtat)] when Xq = by using a population of N particles. 

1. Set n = 0; n is the number of steps. 

2. For i from 1 to N set 

(a) (wH\^\4 ) ) = (0Ax). 
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3. Set 6 = and %, = {(w^ ,x^)} i=lj ... jN . 

4. While the set 9\ n of still running particles at step n is nonempty do: 

(a) Set 9Vfi =0- 

(b) Do #lH n times the following operations: 

(i) Pick a still running particle of index j at random according to 
a family of discrete probability distribution 



Pj ~ T w {k) ' 

/—jk index of particles in dKn n 

where Wn is the weight of the particle after n iterations. 

(ii) The particle is moved in time and space according to the exit 
time and position from a time-space parallelepiped that con- 
tains (tn,Xn ). Its new position is denoted (i^2i> x n+i) and its 
associated weight w^U. 

(iii) If t^L = T or if ajjjjL belongs to an absorbing boundary, then 

( w n+i^n+i^ x n+l) * s a dded to the set & of stopped particles. 
Otherwise, it is added to SHn+x- 

(c) Increment n by 1. 

5. Return 



Ei=iw (l} i= i 

when 6 = {(w«,#,x«)}i =1 ...,iv'- 



(0 x (ty 



As we need to keep track of the positions of all the particles, this algorithm 
is memory consuming. On the other hand, it avoids the multiplication of the 
weights. In addition, this algorithm can be modified in the following way: 
instead of using #5K n particles at step n, it is possible to use N particles, and 
in this case, one has to keep track of the number of still running particles 
and to multiply the weights by the proportion of still running particles. 
The algorithm stops when the proportion of still running particles is smaller 
than a given threshold. This approach can be used, for example, for long 
time simulation, or to estimate rare events, as, for example, in [6, 9, 10, 23]. 

4.6. Estimation of the number of steps. Let us consider now the estima- 
tion of the number of steps. In order to do this we will use the techniques 
employed in [26, 28, 29]. 

In Algorithm 2, we have constructed the Markov chain (Ak)k>o which is 
absorbed when reaching I\ = n {T} x D. 
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For a function a on fl, we set 

Pu{t,x) =E n [u(Ai)|A = (t,x)] and A = Pu(t,x) - u(t,x). 
The operator A is the generator of a Markov chain. 

Lemma 1. If T< +00 and 

E n [e 1 \{9 0) E ) = {t,x)]-t> 1 , 

then 

T-t 



E>|(0 o ,3o) = (*,x)]<l + 

7 

Proof. Consider the problem 

J Av(t, x) = —g(t, x), on Q, 



J \ 1 — / j \ j / 7 w 1 

\u(t,x)=0, on [0,T] x T, 



whose solution is 

'v-l 



u(t,x) = E n 



£</(A fc ) 



,k=0 

We remark that if u and g are well chosen this equality gives a good estimate 
of W[v\. 

Let V(t, x) be the function V(t,x) = (T — i)l(t )3; ) e Q- For in Q, we 

have 

x) = i n [y(^i, Si)| (0 O , H ) = (i, x)] - (T - t) < -7. 
Hence T — t> E n [^£~Q 7) (#o, ^o) = (t, x)] and the result follows easily. □ 

Lemma 2. With the previous notation, for every L > fixed, we have 
supP> > L|(0 O , H ) = (t, x)] < (1 + T - t) exp(-c 7 L/(l + T - t)), 

where c is a constant depending on 7; more precisely c converges to 1 as 7 
decreases to 0. 

Proof. The proof follows from the one of Theorem 7.2 in [28]. □ 

Lemma 3. If T = +00, Q is bounded and 

E n [|x + Hi + c| 2 ] > 7 >0, 
where c is such that min^gg \x + c| > C > 0. Then 

with B > max{7, sup^g \x + c|}. 
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Proof. Let us proceed as in [26]. Choose a vector c such that min^gg \x + 
c\ > C > 0, and set 



V(t,x) 
Thus for B 2 > 7, 



5 2 -|x + c| 2 , if (i,x) el + x g, 

0, otherwise. 



AV(t, x) < \x + c| 2 - E n [|a; + Hi + c\ 2 \(0 o , H ) = (i, x)] < £ 2 - 7 
and the result follows. □ 



5. Numerical examples. We present in this section some numerical ex- 
amples in order to test our algorithm. 

5.1. Speeding up the random walk on squares algorithm. In [28] (see also 
[29]), Milstein and Tretyakov propose a method to simulate Brownian mo- 
tions and solutions of SDEs by using the first exit time and position from a 
hyper-cube or a time-space parallelepiped with cubic space basis. A similar 
method has been previously proposed by Faure in his Ph.D. thesis [11]. This 
method is a variation of the random walk on spheres method. Some authors 
already used random walk on squares and rectangles by using the explicit 
expression of the Green function but without simulating the exit time (see, 
e.g., [35]). One of the main features of our approach is the simulation of 
the couple of nonindependent random variables (exit time, exit position) 
by means of real valued random variables. We have explained in [8] how to 
extend this approach to rectangles and the starting point everywhere in the 
rectangle. This approach is still using only one-dimensional distributions. 
However, by using symmetry properties, we can notice that it is simpler to 
deal with squares centered on the current position of the particle than with 
a rectangle. 

Nevertheless, the computation may be time consuming. We are looking 
now to speed up the computations by using a simple density for the exit 
position. 

Let us consider here the d-dimensional hypercube C = [— 1,1] , and a 
fixed time T > (possibly T = +00). Let B be a ai-dimensional Brownian 
motion. We set t b = inf{t > 0\B t £ C}. Let W be a one-dimensional Brow- 
nian motion. We set = inf{t > 0\W t $ [-1,1]}, R(t) = Pobf^ < t], 
r the density of , S(t, y) = F [W t < y\t < r^] and s(t, y) = 8 y S{t, y) 
the density of Wt given {t < t^_ 1 ^ } . 

Let us note that we can easily switch from C = [—1, l] d to any hypercube 
[— L, L] d after a scaling argument in space and time. Thus, from a numerical 
point of view, we need only to implement the required functions r, s, R 
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and S on [—1, 1]. Analytical expressions for these distribution functions are 
easily deduced from the series presented in the Appendix. 

To simulate the exit time and position from [0,T] x C, we proceed in the 
following steps: 

• Compute the probability (3 = 1 - (1 - R(T)) d that t b < T. 

• With probability /3, decide if {t b < T} happens or not. 

• If {t b < T} happens: 

- For a realization U of a uniform random variable U on [0, 1), set 

t b = R~ i (i-(i-\jpy/ d ), 

which is a realization of t b given {t b < T}. 

- Choose with probability l/2d an exit side (J,e), and set £j = e. 

- For each i = 1, . . . ,d, i 7^ J, set Xi = V^h, where the Uj's are d — 1 
independent realizations of uniform random variables on [0,1). With 
probability 1/2, set & = x% — 1 and with probability 1/2, set & = 1 — %j. 

- Compute the weight 

1 TT s(r B ,&) 



It! 



n 



i=l,...,d,i=t=J 



If {t b > T} happens, then: 

- Set T B = T. 

- For i = 1, . . . ,d, set Xi = where the Uj's are d — 1 independent 
realizations of uniform random variables on [0, 1). With probability 1/2, 
set £i = Xi — 1 an d with probability 1/2, set £j = 1 — Xi- 

- Compute the weight 

1 TT s ( T ^i) 



w 



n 



(t b ,£i, . . . , £d) represent the first exit time and position from [0, T] x C, and 
w is the associated weight. 

For the random walk on squares we can also use the idea proposed in [28] 
and in [11]. This leads to the following algorithm: 

• Compute the probability f3 = 1 - (1 - R(T)) d that t b < T. 

• With probability /3, decide if {t b < T} happens or not. 

• If {t b < T} happens: 

- For a realization U of a uniform random variable U on [0, 1), set 

r B = R- i (i-(i-\jpy/ d ), 

which is a realization of t b given {t b < T}. 

- Choose with probability l/2d an exit side (J,e), and set £j = e. 
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Table 1 



Speeding up the 


random walk 


on squares: experiments 


with 1,000,000 samples 


are used 


Method 


T 


Mean of w 


Variance of w 


Time (s) 


Walk on squares 


0.1 






94 


Imp. sampling 


0.1 


1.0005 


0.28 


3.2 


Walk on squares 


0.2 






82 


Imp. sampling 


0.1 


0.9997 


0.014 


1.8 


Walk on squares 


0.5 






10 


Imp. sampling 


0.5 


0.9999 


0.021 


1.2 


Walk on squares 


1.0 






10 


Imp. sampling 


1.0 


0.9994 


0.017 


1 


Walk on squares 


+oo 






10 


Imp. sampling 


+oo 


0.9998 


0.015 


0.98 



- For each i = 1, . . . ,d, J, draw £j according to the distribution of B % _ B 
given t b ' > T B , where t b% = mi{t > Q\B l £ [-1, 1]}. 

• If {t b > T} happens, then: 

- Set t b = T. 

- For i = l,...,d, draw £j according to the distribution of B l _ B given 

T Bl >T B . 

In both cases, we use tabulated values for R and R~ 1 . In order to simulate 
B\ given t b% > t, we use the rejection method proposed by Faure in [11] for 
t E [0.25,2]. Otherwise, we draw B\ by using the fact that it is equal to 
S" 1 ^,!/) for some random variable U with uniform distribution on [0,1). 
This is the method proposed by Milstein and Tretyakov in [28]. For t > 2, 
the latter method is more efficient than the previous one. For t < 0.2, the 
rejection method may give wrong results. For t close to 0.2, the rejection 
method can be up to 6 times faster than the inversion method, while for t 
close to 2, they are comparable in the computation time. 

If the Brownian motion reaches the side labeled by (1,-1) first at time 
T B , then in order to simulate B\ for % = 2, . . . , d we use a random variable 
with density (j)(x) = 1 + x if x G (—1,0] and 4>(x) = 1 — x if x £ [0, —1). In 
this case, the weights w are close to 1 as we see in Table 1, and the execution 
time is usually divided by 10. For T = 0.1, the variance of w is too high and 
leads to some instabilities. In this case, it is preferable to simulate the exact 
distributions of B? given {T < t b }. 

5.2. Solving a bi-harmonic problem. To test the validity of our approach 
with respect to other algorithms, we consider first an example borrowed in 
[28] (see also [29], page 332). Let D = [— 1, l] 2 , and consider the bi-harmonic 
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equation 

(\A 2 u{x) = l, x£D, 

(18) I u(x) = (f>(x), on dD, 

[ \Au(x) = ip(x), ondD, 

with 

(19) <Kxi3±1) = l+A i 0(± i )X2) = i^l, 

(20) ^ 1)± i ) = i±^i, ^ (± i )X2 ) = i±M. 

After setting f(x) = ^Au(x), (18) may be transformed into the system 

{^Av(x) = 1 on D, with u(x) = tp(x) on <9D, 

2- Au(x) — v(x) = on D, with u(x) = (j)(x) on 5.D, 
whose exact solution is 

ry^ _|_ ry"4 ™2 | ™2 

u(x) = — -, v(x) = — -. 

By Ito's formula, it is easy to show that 

u(x) = E[(f>(x + B t b)} - E[t b ^(x + B t b)\ + 2-E[(t b ) 2 ], 

u(x) =E[V'(x + B T s)] -E[r B ], 

where B is a two-dimensional Brownian motion, and t b is, as above, its first 
exit time from D. 

Here, in contrast with the values presented in [28], we only need to use 
one square, since we are not forced to start from its center. We compare 
the results given by our algorithm (first lines) with the one given by the 
random walk on rectangles (second line) . Each side is chosen uniformly with 
probability 1 /4. The time is drawn by using an exponential random variable 
of parameter 1/(1 — exi) if (i,e) is the exit side. The position is drawn 
uniformly on the exit side. This strategy corresponds in some sense to a 
"naive" and simple way to choose the exit time and position. 

As we evaluate quantities of the form E[/(r s , B t b)], we report the quan- 
tities n n ± 2cr n /- v /n, where fi n is the empirical mean of f{r B ,B T B) with n 
samples, and a n is the corresponding empirical standard deviation. The in- 
terval [n n — 2a n / \fn, ii n + 2a n / -y/ra] represents the 95.5% confidence interval 
for K[f(r B , B t b)]. The estimations u(x) and v{x) of u and v for three points 
are given in Table 2. 

Although a small numerical bias seems to appear, our algorithm provides 
results comparable with the random walk on rectangles method. The execu- 
tion time is much smaller than the one given by this method (also the one 
given by the random walk on squares, for which the simulation of one step 
takes less time, but where more steps are needed). 
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.(a) (b). 

•(c) -(e) 
D .(d) 



S 0.25 



Fig. 2. A simple domain D. 



5.3. Estimation of rare events: Computing hitting probabilities. Let us 
consider the following problem: what is the probability p{x) that starting 
from a point x in a domain D a Brownian motion reaches a part S of the 
boundary dD? It is well known that p is the solution of the Dirichlet problem 

(21) ±Ap(x) = 0onD and p{x) = j £ 

We illustrate our method on the simple two-dimensional domain D drawn 
in Figure 2 and we compute the value of p at the five points marked, respec- 
tively, by (a), (b), (c), (d) and (e) on Figure 2. 



Table 2 

Solution of the bi-harmonic equation: the first line of each row contains the results for 
our algorithm, and the second line contains the results for the random walk on rectangles 



X 


n 


u(x) 


u(x) 


v(x) 


v(x) 


Time (s) 


(0.3,0.5) 


10 4 


0.00588 


0.0047 ± 0.0037 


0.17000 


0.1638 ±0.0081 


0.03 








0.0064 ± 0.0039 




0.1684 ±0.0081 


3.8 




10 5 




0.0061 ±0.0012 




0.1669 ±0.0026 


0.23 








0.0062 ±0.0012 




0.1679 ±0.0026 


38 




10 6 




0.0059 ± 0.0004 




0.1698 ±0.0008 


2.2 








0.0059 ± 0.0004 




0.1696 ±0.0008 


381 


(0.7,0.8) 


10 4 


0.05414 


0.0480 ±0.0017 


0.56500 


0.5297 ±0.0064 


0.02 








0.0553 ± 0.0020 




0.5707 ±0.0061 


7 




10 5 




0.0526 ± 0.0005 




0.5593 ±0.0019 


0.2 








0.0543 ± 0.0006 




0.5652 ±0.0019 


73 




10 6 




0.0536 ± 0.0002 




0.5654 ±0.0006 


2.5 








0.0542 ± 0.0002 




0.5650 ±0.0006 


726 


(0.9,0.9) 


10 4 


0.10935 


0.1103 ±0.0009 


0.81000 


0.8186 ± 0.0034 


0.01 








0.1109 ±0.0020 




0.8105 ±0.0038 


11 




10 5 




0.1131 ±0.0002 




0.8390 ± 0.0006 


0.2 








0.1095 ±0.0003 




0.8107 ±0.0011 


112 




10 6 




0.1087 ±0.0001 




0.8097 ±0.0003 


2 








0.1093 ±0.0001 




0.8100 ±0.0003 


1100 
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Table 3 

Computation of p(x) at given points of D 



Point Import, sampling Finite element Walk on rect. 



(a) 


3.32 


1(T 6 


3.39 


1(T 6 


0.00 




(b) 


2.31 


1(T 5 


2.23 


1CT 5 


1.00 


10~ 5 


(c) 


1.70 


io- 4 


1.77 


10~ 4 


1.90 


10~ 4 


(d) 


4.43 


10~ 5 


4.64 


10~ 5 


3.00 


10~ 5 


(e) 


2.79 


10" 3 


2.81 


10" 3 


2.36 


10" 3 



To set up our algorithm, we use two rectangles as in Figure 3. The numbers 
marked on each side are the probabilities to reach each one of these sides. 

In order to obtain the simulated exit time we draw an exponential random 
variable with parameter a where a is given by a = l/(y Li/2). The Lj notes 
the length of the rectangle in the direction perpendicular to the boundary 
that the particle hits. 

We perform 100,000 samples; each computation takes around 1 s on our 
computer (a MacBook 12", 2 GHz with a code written in C). The values for 
p are given in Table 3. We perform a comparison with the value given by 
MATLAB/PDEtool where (21) is solved by using a finite element method, 
and with the method of random walk on rectangles [8] which is exact (up 
to the Monte Carlo error), for such a domain. In this case, with a sample of 
size n, the variance of the empirical mean is p{x){\ — p(x))/n. 

We notice that the results given by our method are close to the one 
given by the finite element method. As one can expect, the random walk on 
rectangles (and any other methods that do not rely on importance sampling 
or variance reduction techniques) is not efficient to estimate the values of 
p(x) when they are of the same order as the standard deviation of the 
empirical mean. 

In order to test the validity of our method for the simulation of rare 
events, we use the domain D' as in Figure 4. 



1/10 







1/4 


1/10 


1/4; 

i_ 








1/4 






7/10 


1/10 





Fig. 3. Decomposition of D into rectangles. 
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The numerical results are reported in Table 4. p n is the empirical mean 
with n = 100,000 samples, and S5o(Pn) is the empirical standard deviation 
computed over 50 realizations of p n . We obtain really good results even while 
computing small probabilities of order of magnitude 10~ 10 . 



5.4. Simulation of SDEs: Approximation close to the boundary. Let us 
consider the two-dimensional SDE 

1 \s\xi(xi+ x 2 ) 
1 



X t = x+ [ a(X s )dB s with a(x) 
Jo 



which is driven by a two-dimensional Brownian motion B. The process X 
is killed when it exits from the domain D which is represented in Figure 5. 

In order to simulate X, we use either an Euler scheme with a time step of 
0.0025 or a (possibly modified) random walk on squares. The squares sides 
lengths are smaller than 2L with L = 0.05 (note that the time step of the 
Euler scheme corresponds to 0.05 2 which is close to the average exit time of 
the square [0.1, 0.1] 2 ). As the diffusion moves in a bounded domain, we use 
to deal with the boundary condition and apply the technique proposed in 
[7]: if the distance between the position of the particle and the boundary is 
smaller than 2L, we choose the square such that one of its sides is included 
in the boundary when it is possible to do so. 

Unless the coefficients of the SDE are constant, one needs to simulate 
many couples of exit times and positions from small squares, and the com- 
putational time becomes very large and is not competitive with respect to 
the Euler scheme. In addition, when the random walk on squares is coupled 
with importance sampling, the weights grow quickly (see Section 4.4). 

When the Euler scheme is used, we simply stop the algorithm when the 
particle leaves the domain D. This is a crude way to proceed, and some 
refinements can be done (see, e.g., [15]). Note that the exit time is then 
overestimated. 

The idea is to mix the two methods and to use the Euler scheme inside 
the domain, and a random walk on squares when the particle is close to the 
boundary. We improve thus the simulation as in this case the behavior of the 
particle is taken into account. In addition, it is possible by making a change 



•(a) (b). 

• (c) 
D' .( d) 



(e) \ S (width 0.1) 



Fig. 4. A simple domain D' . 
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Table 4 

Computation of p(x) at given points of D' 



Point 


Pn 




Finite element 


(a) 


1.00- lO -10 


2.3 -lO -11 


1.15- 1(T 10 


(b) 


7.67- 1(T 10 


1.6 -10 -10 


8.13- 10" 10 


(c) 


5.19- 1(T 9 


1.0- 1(T 9 


6.61 • 10 -9 


(d) 


1.31 ■ 1(T 9 


2.8- 1(T 10 


1.73- 10 -9 


(e) 


2.27- 1CT 7 


4.9- 10~ 8 


2.29- 10" 7 



of measure, to increase or to decrease the probability that the particle hits 
the boundary. 

Our aim is here to increase the number of particles which are not killed 
before a given time T. When one side of the square is set on the boundary, 
we use a probability p that the particle reaches the side of the square that 
is opposite to the boundary, and q = (1 — p)/3 for any other side. We have 
thus a "repulsing" effect. 

We use P 1 = {p = 0.7, q = 0.1} and P 2 = {p = 0.91, q = 0.03}. 

In order to avoid the explosion of the variance of the weight, we have 
used a limitation N max for the number of times this procedure is used. The 
variance of the weight for each time this procedure is used is 0.04 for the set 
Pi and 0.34 for the set P 2 . 

All the simulations are done with 100,000 particles. The results are sum- 
marized in Table 5. For T = 1, the proportion of particles still alive is of order 
0.19% (using the Euler scheme without specific treatment on the boundary, 
we get an estimation of 0.33%, yet for a quicker simulation of 7 s). With a 
population Monte Carlo method, we obtain an estimate of 0.17%, using the 
set P\ and a running time of 126 s. We see that our scheme allows one to 
get much more alive particles. 



1 unit 




Fig. 5. Domain D with the label of the sides and the starting point. 
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Table 5 

Simulations of the proportions (in %) of the particles reaching a given part of the 
boundary as well as the surviving particles at time T. We write "unstable" in the column 
for the variance of weights when the mean of the global weights is rather far from 1 



T 


Type 


Side 1 


Side 2 


Side 3 Side 4 


Side 5 


Final time Var. weights 


jV max 


Time 








Test 


with 


set of probabilities Pi on 


the boundary 






1 


Est. 


31.55 


12.39 


4.16 


0.44 


51.27 


0.17 


9.3 


5 


73 




Sim. 


31.91 


12.97 


5.07 


0.65 


48.96 


0.41 








1 


Est. 


30.81 


13.19 


4.06 


0.32 


51.43 


0.17 


29.9 


10 


83 




Sim. 


31.93 


13.08 


5.42 


0.75 


48.10 


0.66 








1 


Est. 


31.05 


13.83 


4.37 


0.42 


50.14 


0.17 


30.0 


20 


93 




Sim. 


32.01 


13.40 


5.57 


0.91 


47.10 


0.96 








1 


Est. 


30.96 


13.54 


4.08 


0.36 


50.84 


0.19 


56.55 


100 


99 




Sim. 


31.78 


13.27 


5.57 


0.98 


46.83 


1.36 














Test 


with 


set of probabilit 


ies P2 on 


the boundary 






1 


Est. 


29.45 


12.11 


3.49 


0.58 


54.19 


0.14 


426 


5 


90 




Sim. 


32.13 


13.07 


5.71 


0.81 


47.61 


0.95 








1 


Est. 


33.76 


11.78 


5.71 


0.37 


48.16 


0.19 


65.5 (unstable) 


10 


117 




Sim. 


32.03 


13.50 


6.70 


1.13 


45.21 


1.48 








1 


Est. 


31.28 


14.19 


3.75 


0.44 


50.10 


0.21 


51.08 (unstable) 


20 


162 




Sim. 


31.18 


13.48 


7.64 


1.62 


42.44 


3.62 








1 


Est. 


29.87 


13.73 


2.83 


0.30 


53.03 


0.23 


312.5 (unstable) 


100 


223 




Sim. 


28.13 


12.21 


7.50 


1.58 


36.36 


14.23 









APPENDIX: HOW TO GET DENSITIES FOR DIFFERENT 

SITUATIONS? 

We present in this section analytical expressions for the density in different 
cases. 

Except for the case of a drifted Brownian motion with Dirichlet boundary 
condition at one endpoint of [—1, 1] and a Neumann boundary condition at 
the other endpoint of [—1,1], we obtain two expressions, one which follows 
from the images method and the other one from the spectral decomposi- 
tion. From a numerical point of view, the spectral decomposition gives rise 
to series that converge very quickly for large times. It is worth using the 
expressions given by the method of images for small times. 

A.l. Brownian motion without drift. We are interested in this section 
in writing down some useful formulas for the calculations. Let us consider 
first the case of the standard one-dimensional Brownian motion starting 
from x 6 [—1, 1] which is killed or reflected when hitting the boundaries — 1 
or 1. We shall write D for Dirichlet condition on the boundary and N for 
Neumann condition, which of course correspond to killing and, respectively, 
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reflection. Furthermore we shall note, for example, PDN(t, %i, x 2 ) the density 
of the Brownian motion on [—1, 1] killed when hitting — 1 and reflected on 
1 more precisely the order in the indices indicates the boundary condition 
in —1 and 1, respectively. 

A. 1.1. Reflected Brownian motion on [—1,1]. Let PNN{t,x\,x 2 ) denote 
the probability density function of a Brownian motion at time t, starting 
from x\ and reflected at —1 and 1. By using the method of images we get 
the following formula for the transition density: 



^ 00 

PNN(t,X 1 ,X 2 ) = -L= V [ e -(^-2+4n)V(2*) +e -( :E1+ a :2 +4n+2)V(2t)]. 
ylut 

v n=— oo 

The spectral representation of this density writes 



PNN(t,xi,x 2 ) = ^ + £e ^ 2 / 8t cosfcx 1 + l)) COS 

n=l ^ ' 

These expressions may be found, for example, in [4]. 



(-2-0=2 + 1). 



A. 1.2. Killed Brownian motion on [—1, 1]. Let PDD(t,x\,x 2 ) denote the 
probability density function of a Brownian motion at time t, starting from 
x\ and killed when it exits from the interval [—1,1]. That is, 

PDD(t,xi,x 2 )dx 2 =F Xl [B t € dx 2 ;t< tdd], 

where tdd — inf{i > 0; Bt ^ [—1, 1]}. Then, by the images' method we have 



1 00 

PDD(t,X 1 ,X 2 ) = -L= V [ e -(-i-2+4n)V(2t) _ e -(. 1+ . 2+ 4n + 2)V(2: 

- ' 2nt 



n=— oo 



For the law of the exit time we get 

^ 00 

Pa* [tdd G dt] = -_ (-l) n (^i + 2n + l) e -^+ 2 « +1 ) 2 /(2*) dt 

v n=—oo 

The spectral representation can be also written and yields 

PDD(t, Xl ,x 2 ) = f; e -« 2 - 2 /8* sinfc^ + 1)) sinfe* 2 + 1) 

n=l ^ ' ^ 

The law of the exit time is given by 
P*! [tdd €dt] = ^ f;(-l) re (2n + i) e -(2-+i) 2 - 2 /8 < cos (L+^\ nxA dt. 

n=0 ^ ' ' 

These expressions may be found, for example, in [4] or in 
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A. 1.3. Mixed boundary conditions for the Brownian motion on [—1,1]. 
We give here explicit solutions for the Brownian motion killed on —1 and 
reflected on 1. Let pDN(t, xi, x 2 ) denote the probability density function of 
a Brownian motion at time t, starting from x\ and killed when it hits — 1 
and reflected on 1. Then, by the images' method, one gets 



00 

p DN (t,X 1 ,X 2 ) = -L= V (_l)»[ e -(*i-*a-Mn) a /(2t)_ e -(*i+* a -M«+2)V(a 

v 2nt 



n=—oo 



Let us denote also by tdn the killing time for the Brownian motion on 
[—1,1) killed on —1 and reflected on 1. Hence 

00 

¥ xl [r DN edt] = ^L= V (_i r(a;i+ 4n + l)e-^ +4 " +1 ) 2 /( 2 *)dt. 
V 2vr^ n= _ oc 

The spectral representation can be also written and yields 

PDN (t, Xl ,X 2 ) = jr e -(*H-D^/32t gin f^±l)l {xi + 1)) 
n=0 



/(2n + l)vr, 
x sin I (x 2 + 1, 

Then we get from the spectral representation the law of this exit time, 
F X1 [r DN Edt] = l f>n + l)e-( 2 " +1 ) 2 - 2 / 32t sin( (2n + 1)7r (x 1 + 1)) dt. 

n=0 ^ ' 

The dual situation (reflection on —1 and absorption on 1) can be obtained 
easily by the transformation 

PND(t,Xi,X 2 ) =PDN(t,-X 1 ,X 2 ). 

These expressions may be found, for example, in [4]. 

A. 2. Brownian motion with drift fi. As in the previous part of the 
Appendix we consider here the case of the Brownian motion with drift on 
the interval [—1,1] which is killed or reflected on —1 and 1. If we note by 
p..'^(t,xi,x 2 ) the law of the process with drift fi and living on [— L,L] and 
p~(t, x\, x 2 ) the corresponding law on [—1,1], then by the properties of the 
Brownian motion we have 



where the dots in the indices can take the value D for a Dirichlet condition 
or N for a Neumann condition, as previously noted. 
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A. 2.1. Brownian motion with drift \i reflected on [—1,1]. We keep the 
same notation as before. The use of the images' method gives the following 
representation of the density: 



2fie 2flX2 

= 2/i p— 2fi 



+ 



1 



'2irt 



E 



e 4/m e -(xi-x 2 +M+4n) 2 /(2t) 



+ 



oo 

^= E 

/ 27rt_ . 



e -2tiX! e -fi(4n+2) e -(x 1 +x 2 -lit+4n+2) 2 /(2t) 



n=—oo 
oo 



^2 J- e^ 4 " +2 )erfc 



2Cl + x 2 + A 4 * + ^ n + ^ 



This formula can be obtained also from the results in Veestraeten [37]. 
By the spectral method (see, e.g., [24]), we have, after some calculations, 



2/ie 2 ^' 2 

g2fi g— 2/1 

li(x2-xi)-fi 2 /2t 



-n 2 w 2 /8t 



E 



^p? + n 2 7r 2 /4 



7rn / nn 



( 7rn 1 
+ //sinl — (xi + 1) 



/ 7T 



j cosj^— (x 2 + l)j +/xsin^— (x 2 + r 



A. 2. 2. Brownian motion with drift fi on [—1,1] killed at the boundary. 
We keep the same notation as before. By using classical properties of the 
Brownian motion and the results from Milstein and Tretyakov [28] we have 
the following transformation: 

p fl DD (t,x 1 ,x 2 ) = e^ X2 - x ^ 2t / 2 PDD (t,x 1 ,x 2 ). 

Then, by the images' method, 



Li(x2—xi)—[i 2 t/2_ 



'2%t 



E 



J e -(xi-x 2 +4n) 2 /(2t) _ e -(xi+x 2 +4n+2) 2 /(2t)j_ 
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We write down both distribution and density for the exit time. The distri- 
bution writes 



OO r- , 

nK D <t] = l-- £ ^ erfc 



' x\ + fit + An — 1 



erfc ( 



2t 



x\ + fit + 4n + 1 



2t 



1 OO 

+ _ e -(2Atxi+/x(4n+2)) 



. f Xi - [it + 4n + 1\ 

errc — 

V V2i J 



erfc 



x\ — \xt + 4ra + 3 



2/ 



while for the density we obtain 

-fix 1 -fi 2 t/2 00 



Sdt] = — -=^- V [e-^(xi+4n + l)e^ (:Cl+4n+1) ' 



/(2t) 



_ e M( Xl + 4 n _l) e -(*i+4n-i)V( 2 *)] 
The spectral representation can be also written and yields 

p fl DD (t,x 1 ,x 2 )=e^- x ^ 2t / 2 



( nn 



n=l 



The distribution of the exit time is given by 
PxJrg D <t] 

-fj,xi—fi?t/2'S~ y (—fi ■ — - 



1-e" 



^E(e^-(-l) 



n=l 



4/x 2 + n 2 7r 2 



-n 2 7r 2 /8t 



, n7r 

x sm( -^-{xi + 1, 



n=l 

oo 



^2 _|_ n 2^2 



( e -M + e M)y-(_i r 2(2n + l)7r 
1 + ^ J 4u 2 + (2n+l) 



n=0 



x e 



4/x 2 + (2n + l) 2 vr 2 

(2n+l)^/8f CQ / (2n+l)7T 



-XI 
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and 



Bl b% D €dt] = e-'™-'> l ' 2 



oo 
n=l 4 



x sin^^(xi + 1) j dt. 
In a more detailed expression we can write this on the form 

oo 

x J2(-l) n —e- n2 * 2 / 2t S m(mr Xl ) 

n=l 



x V( 1 - ) n ( 2n + 1 ) 7r c -(2n+l) 2 7r 2 /i 
n=0 



/(2n+l)vr \ 
x cos I xi J dt. 

These expressions may be found, for example, in [4] or in [28]. 

A. 2. 3. Mixed boundary condition for the Brownian motion on [—1, 1] with 
drift The aim is to express some explicit solutions for the Brownian 
motion killed on — 1 and reflected on 1 . We solve now the following eigenvalue 
problem: 

\<ff'(x\) + (J,ip'(xi) = Xif(xi), 
¥>(-!) = 0, 
¥>'(!) = 0- 

We can remark first that if <p\ is an eigenfunction for the eigenvalue A for 
the preceding PDE, then A is negative. 

We associate with this problem the corresponding second degree equation 
and note A = /i 2 + 2A. After a detailed calculus with respect to the sign 
of A we can express the countable set of eigenfunctions and eigenvalues 
with respect to the possible values of \x. There are three different situations, 
expressed in Table 6 (see, e.g., [33]). The density j3dat(£, 3:1,22) is obtained by 
using the spectral expansion PDN(t,Xi,x 2 ) = 52k>o ex P^kt<P\ k (xi)ip\ k (x 2 ), 
where • • • < A2 < Ai < Aq- The density qDN(t,xi) of the exit time is also 
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Table 6 

Eigenvalues and eigenfunctions for the Dirichlet/ Neumann problem 
with a constant transport term fi 



/x A ip\ 

c -i"i 

N /2(l-(co S 2(2- s /3^2Z5I)) / (2^)) 



V~^ 2 -2A 



tan(2 v /-^i 2 -2A) = 

^ = 2 4 

A<-± - e ' X ^ -sin(,/i+2A(l + a;i)) 

8 y/2\ sin(2^/l/4+2A) v V 4 

tan(2 ^/(|+2A)) = 2 ^/(|+2A) 

\J2 cosh 2 (2-^/ n 2 +2A)/f-l 



^-^(xi + 1) 



M >± *>-%-, ^ - : -M /.- , i °'h i l)) 



tanh(2VM 2 + 2A) = Vm ' +2A ^f* 1 sin( \/-^ 2 - 2A>1 + 1)) 

M \Z2(l-co S 2(2y- M 2_2A)/(2n)) 

A<-4 . 
tan(2 x /-^ 2 -2A) = V^ 2 ~ 2 * 



expressed by 



ei [tew £ = - y^\ k e Xkt (f)\ k (xi) / <t>\ k (x 2 ) dx 2 . 

k>0 
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